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This  work  is  part  of  an  ongoing  research  program.  New  developments,  refinements,  and  correc¬ 
tions  evolve  continuously,  and  many  specifics  in  the  report  may  have  been  superseded  by  the  time 
it  is  read.  Readers  are  encouraged  to  communicate  with  the  authors  to  obtain  more  information 
on  latest  developments. 


NUMERICAL  SOLUTIONS  FOR  A 
RIGID-ICE  MODEL  OF  SECONDARY 
FROST  HEAVE 

Kevin  O’Neill  and  Robert  D.  Miller 


INTRODUCTION 

Primary  frost  heave  occurs  when  ice  grows  on  the  edge  of  some  soil.  The  ice  is  fed  by  a  liquid 
flow  through  the  contiguous  soil,  but  the  ice  does  not  significantly  penetrate  the  soil.  In  the  more 
common  instance  called  secondary  heave,  heave  proceeds  when  ice  has  penetrated  the  soil.  A  rigid- 
ice  model  of  secondary  frost  heave  in  air-free,  solute-free,  colloid-free  soils  has  been  described  both 
qualitatively  (Miller  1977)  and  quantitatively  (Miller  1978),  but  the  interacting  equations  constituting 
the  quantitative  model  have  been  solved  only  for  very  simple  quasi-steady  states  (Miller  and  Koslow 
1980).  In  this  paper  we  report  test  case  solutions  of  the  equations,  simulating  the  course  of  frost 
heave  in  a  hypothetical  soil  column  under  a  specified  load. 

The  physical  theory  relied  upon  here  is  that  in  the  references  above.  Its  salient  features  include 
nonisothermal  freezing,  such  that  temperature,  pressure,  and  ice  content  are  all  connected  by  a 
single  thermodynamic  relation,  and  each  of  those  quantities  varies  continuously  across  the  freezing 
zone.  In  addition,  unfrozen  moisture  flows  across  this  freezing  zone,  which  is  of  finite  thickness 
("frozen  fringe”).  This  flow  feeds  ice  lens  growth,  which  is  accompanied  by  large,  continuous 
changes  in  pressure  and  hydraulic  conductivity  across  the  fringe.  An  ice  lens  grows  because  1 )  pore 
ice  is  "extruded"  upward  by  thermally  induced  regelation  and  2)  ice  forms  from  water  moving  up¬ 
ward  through  unfrozen  films  between  the  pore  ice  and  the  soil  particles  immediately  beneath  the 
lens.  Once  freezing  begins,  pore  ice  forms  by  accretion  on  pre-existing  ice.  Thus,  all  portions  of 
ice  are  rigidly  interconnected  and  move  as  a  rigid  body.  The  theory  develops  criteria  in  terms  of 
thermodynamic  variables  for  the  sequential  locations  of  lenses  as  freezing  progresses. 


GOVERNING  EQUATIONS 

The  equations  used  here  to  simulate  freezing  and  heaving  processes  come  from  the  general  laws 
of  mass  and  energy  conservation,  together  with  thermodynamic  and  other  specialized  equations 
appropriate  for  saturated  freezing  soils.  Of  central  importance  is  the  expression  for  the  pressure 
jump  due  to  curvature  in  the  tiquid-water/ice  interface: 


U)~U  = 
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where 


u  =  liquid  pressure  (N/m2) 
u j  =  ice  pressure  (N/m2) 

Wiw  =  ice/liquid-water  interfacial  tension  (0.0331  N/m;  Koopmans  and  Miller  1 966) . 

The  parameter  i//(m-1 )  may  be  interpreted  as  a  mean  curvature  of  the  ice/liquid  interface.  As 
some  soil  becomes  progressively  more  frozen,  ice  penetrates  smaller  and  smaller  pores,  and  the 
average  curvature  of  the  water  phase  interface  increases. 

An  equation  identical  to  eq  1  holds  for  the  pressure  jump  across  the  air/liquid-water  interface 
in  unfrozen,  unsaturated  soil,  with  air  pressure  substituted  for  ux  and  air/water  interfacial  tension 
substituted  for  wiw.  Pursuing  this  similarity,  Koopmans  and  Miller  have  shown  that  for  a  given 
soil  a  “freezing  characteristic  curve”  may  be  obtained,  relating  volumetric  ice  content  6 to  \p.  This 
curve,  in  turn,  is  related  to  the  soil  moisture  characteristic  curve  simply  by  a  scaling  constant.  In 
addition  the  parameter  \j/  may  be  expressed  using  a  Clapeyron  relation: 

\p  =  Au  +  BT  (2) 

where  T  denotes  temperature  (°C),  and  the  constants  A  and  B  are  -2.51  m/N  and  -3.40x  107 
(m°C)_1 ,  respectively.  Thus 

fl,  =  fli(0)  =Q.,(Au  +  BT).  (3) 

One  can  use  eq  3  to  express  the  differential  of  8X  as 


To  construct  a  governing  set  of  differential  equations,  one  may  begin  by  considering  conserva¬ 
tion  of  mass  over  a  small  volume  element  of  soil,  obtaining 

^  [p0  +  p ,0,]  +  ^  [pv  +  Pji/,]  =0  (5) 

where 

0  =  volumetric  liquid  content 
p  =  density  of  liquid  water  (1000  kg/m3) 
v  =  liquid  volume  flux  (m/s) 

Pi  =  ice  density  (917  kg/m3) 

Vj  =  ice  flux  (m/s) 
t  =  time(s) 

x  =  space  coordinate  (m),  positive  downward. 

The  ice  flux  is 

•'l  *  (6) 

where  the  ice  velocity  Vt  (m/s)  is  variable  in  time  but  constant  in  space,  in  keeping  with  the  rigid- 
ice  assumption.  We  assume  that  liquid  flux  may  be  expressed  in  both  freezing  and  unfrozen  zones 
using  a  Darcy-type  law: 
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where  p  is  the  acceleration  of  gravity  (9.8  m/s2)  and  A(0j)  is  the  hydraulic  conductivity  (m/s). 
Noting  that 


0j  =  0o-0 


(8) 


where  0O  is  the  soil  porosity,  and  incorporating  eqs  6  and  7,  one  can  rewrite  eq  5  as 


be, 

<*■-'>  aT 


=  o. 


(9) 


For  the  most  part,  it  seemed  best  to  use  eqs  3  and  4  to  translate  all  expressions  containing  6 
and  0j  explicitly  into  expressions  in  u  and  T.  One  can  then  express  most  physically  reasonable 
boundary  conditions  relatively  simply.  Also,  the  same  set  of  governing  equations  in  u  and  T  may 
be  applied,  with  appropriately  specified  parameters,  in  both  the  unfrozen  and  the  frozen  zone. 
To  this  end,  bd-Jbt  in  eq  9  may  be  expressed  using  eq  4,  yielding 
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where  A p  denotes  (pj-p). 

Energy  conservation  within  a  small  representative  volume  of  soil  provides  another  equation: 
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where 


cn  =  heat  capacity  of  component  n  (J/kg) 

Kh  =  soil  thermal  conductivity  (J/m  s°C) 

L  =  latent  heat  of  fusion  of  water  (3.35x  10s  J/kg). 

The  terms  bO  Jbt  +  V\b8  Jbx  are  the  rate  at  which  ice  volume  forms  per  unit  volume  of  soil. 

The  formation  rate  is  not  simply  bO  Jbt,  because  ice  may  accumulate  within  a  portion  of  the  soil 
due  to  regelation  heave  processes,  that  is,  due  to  migration  as  opposed  to  net  phase  change.  The 
rate  of  phase  change  is  equal  to  the  difference  between  the  ice  accumulation  rate  and  the  in-migra¬ 
tion  rate,  that  is,  to  the  sum  of  d0j/9f  and  Vjbd-Jbx.  The  particular  form  of  these  terms  arises  from 
a  systematic  derivation  of  the  equation,  not  merely  from  these  remarks. 


Using  eq  4  to  reexpress  d 0,/ 
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SOLUTION  PROCEDURES 

If  the  remaining  terms  in  eqs  1 0  and  1 2  containing  V\  are  reexpressed  using  eqs  3  and  4,  then 
the  former  set  becomes  two  equations  in  the  unknowns  u  and  T.  In  practice,  the  Vi  terms  are  in¬ 
cluded  in  the  overall  u  and  T  solution  iteratively.  The  influence  of  these  terms  is  relatively  mild, 
so  that  for  any  current  set  of  u  and  T  values,  0j  expressions  in  the  equations  can  be  approximated 
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for  the  immediate  future.  This  approximation  is  then  improved,  as  the  equations  are  re-solved  at 
the  same  point  in  time  for  more  consistent  u  and  T  values. 

A  finite  element  scheme  is  used  for  the  solution  in  space,  with  simpler  finite  differences  in  time. 
This  is  frequently  done  in  solving  transport  equations  to  treat  higher  derivatives  in  space  and  strongly 
space-dependent  parameters  accurately.  This  method  has  been  found  to  be  flexible  and  efficient 
for  simpler  freezing  problems  (Lynch  and  O’Neill  1981).  In  the  finite  element  system,  dependent 
variables  are  expressed  as  finite  sums  of  predetermined,  space-dependent  interpolation  ("basis") 
functions,  each  multiplied  by  its  coefficient,  which  is  time  dependent  in  this  case: 

N  N 

T=Z  /■jMH'jM;  u  =  2  Ux(t)W;(x) .  (13) 

i  i  i=l 

In  this  study  the  basis  functions  H-'j  provide  linear  interpolation  of  the  unknowns  between  mesh 
points,  and  the  coefficients  Tf  and  U-t  correspond  to  dependent  variable  values  at  the  mesh  points. 
The  problem  is  solved  by  evaluating  all  unknown  coefficients  in  eq  13  through  time. 

Substituting  eq  13  into  eqs  10  and  1 2  and  using  well-established  Galerkin  finite  element  pro¬ 
cedures  (e.g.  Pinder  and  Gray  1977)  results  in  two  coupled  governing  algebraic  (matrix)  equations, 
each  of  the  general  form 

I CJ  M+  10)  »!/}  =|/?}  (14) 


where 


)l/j  =  vector  of  unknowns,  composed  of  both  temperature  and  pressure  coefficients 
=  its  time  derivative 

|/?|=  vector  containing  specified  quantities  and  the  V,  terms. 

At  this  point  the  problem  has  been  transformed  into  a  matrix  ordinary  differential  equation  in  time. 
Time  derivatives  of  the  coefficients  in  {l^fare  expressed 


it+1  -  yf 


(15) 


where  A f  is  the  time  step  size  and  k'-1  denotes  the /th  coefficient  evaluated  at  the  fcth  time  level. 
The  time  steps  may  vary  in  size,  ranging  in  this  study  from  seconds  as  freezing  begins  to  hours  as 
a  steady  state  is  approached. 

Using  eq  15  in  eq  14  results  in  two  matrix  equations  of  the  form 


(£]  |lA+,(  =]/?}-  (£)  . 


(16) 


In  this  equation  the  matrices  [£|  and  [£|  are  time  dependent  because  they  contain  quantities  de¬ 
rived  from  the  dependenrvariables  (e.g.  ddjd\p).  Optimal  stability  was  obtained  with  a  "fully 
implicit"  formulation,  in  which  the  matrices  are  evaluated  at  the  £+1  time  level,  and  eq  15  repre¬ 
sents  a  backwards  difference  in  time.  Thus,  in  the  course  of  a  simulation  the  matrices  are  updated 
iteratively.  The  boundary  conditions  in  terms  of  the  coefficients  are  also  incorporated  in  eq  16. 

In  general,  the  solution  cycles  may  be  characterized  as  follows:  Both  temperature  and  pressure 
coefficients  in  the  l/'j5  are  known,  either  from  previous  analysis  or  from  the  initial  conditions.  On 
the  basis  of  these  values,  dO-Jdip,  fc(0|),  etc.  are  estimated,  from  which  the  matrices  and  /?|  at  the 
next  time  level  are  specified.  The  system  is  then  solved  for  the  values,  i.e.  for  the  dependent 
variable  values  at  the  £+1  time  level.  These  values  are  then  used  to  re-estimate  other  dependent 
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quantities  at  1 .  The  system  is  then  solved  again  for  a  new,  presumably  more  consistent  set  of 
1  values.  The  process  is  repeated  until  the  degree  of  convergence  is  deemed  sufficient  for  one 
to  proceed  to  the  next  point  in  time. 

At  any  point  in  time,  a  formula  for  l/|  is  obtained  in  terms  of  the  current  values  of  u  and  T  by 
taking  a  mass  balance  at  the  bottom  of  the  lowest  lens.  Throughout  the  lens  (including  its  base) 
ice  is  moving  upwards  with  a  mass  flux  of  Just  below  the  lens,  soil  ice  is  rigidly  connected 
to  the  lens  ice  and  therefore  has  the  same  velocity  Vt.  However,  below  the  lens,  ice  composes  only 
a  fraction  of  the  medium  cross  section  equal  to  6 ^  thus,  the  ice  flux  there  is  p-,  l/, j.  A  liquid  mass 
flux  pv  also  exists  in  the  soil  below  the  lens.  Equating  the  mass  flux  entering  the  lens  boundary 
from  below  with  that  leaving  from  above  yields 

pil/,  =  pil/,0i  +  pt'.  (17) 

When  eq  7  is  substituted  into  eq  1  7,  one  obtains  an  expression  for  V{: 


When  k  and  0{  arc  evaluated  from  expressions  in  u  and  T,  Vf  may  be  evaluated  for  use  in  the  gov¬ 
erning  equations  and  in  other  calculations. 


THE  COURSE  OF  A  SIMULATION 

So  far,  simulations  using  the  theory  in  this  paper  have  been  performed  for  a  column  15.3  cm 
long.  The  warm  end  was  held  at  1°C  and  atmospheric  pressure  (u  -  0);  initially  the  entire  column 
was  considered  to  be  uniformly  at  1°C  and  zero  gage  pressure.  The  temperature  on  the  freezing 
side  soil  surface  was  specified  through  all  time.  The  boundary  conditions  at  the  freezing  end  varied 
according  to  the  stage  of  the  simulation,  as  described  below. 

In  the  first  stage  the  freezing  side  surface  temperature  is  lowered  at  a  constant  rate  from  0°C. 
With  a  specified  overburden  pressure  P ,  heave  is  temporarily  restrained,  and  the  ice  velocity  K,  is 
therefore  zero.  The  liquid  flux  at  the  freezing  end  is  also  zero  so  that  one  can  specify  a  pressure 
boundary  condtion  there  using  eq  7: 


(dulbx)x=Q  =  pg  .  (19) 

During  this  stage  the  ice  content  increases,  driving  liquid  out  of  the  column,  and  the  liquid  and  ice 
pressures  both  increase. 

The  second  stage  of  the  simulation  begins  when  the  ice  pressure  at  the  freezing  surface  (as  cal¬ 
culated  from  eq  1 )  reaches  P,  so  that  a  lens  there  can  support  the  overburden.  Heave  and  lens  for¬ 
mation  begin  and  V (  is  nonzero.  The  analysis  concentrates  on  the  freezing  zone  below  the  lowest 
lens,  where  the  crucial  heat  and  mass  transfer  activities  take  place.  The  mesh  point  separations  are 
as  little  as  2.5x  10-5  m  immediately  below  the  lens,  and  they  increase  geometrically  to  about  3  cm 
at  the  unfrozen  end.  The  freezing  activity  above  the  lowest  lens  is  assumed  to  be  slight,  so  that 
the  temperature  distribution  there  may  be  considered  linear.  This  assumption  was  borne  out  by 
results  in  the  mostly  frozen  zone  before  heave  began.  (This  assumption  was  enlisted  only  as  a  com¬ 
putational  convenience  for  the  current  examples  and  need  not  be  employed  in  more  sophisticated 
simulations  based  on  the  same  theory.) 

One  boundary  condition  for  the  soil  below  the  lowest  lens  is  obtained  from  a  heat  balance 
across  the  lens  boundary: 

*hf  (rLr±)=  (|~)x=b  -  Lpv(*=t>)  (20) 
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where 


(x  =  b)  =  location  of  the  base  of  the  lowest  lens 

8  =  length  of  the  frozen  material  between  the  soil  surface  and  x  =  b 
T%  (t)  =  specified  surface  temperature 
Tb  ( t )  =  temperature  at  x  =  b 

Kh f  =  thermal  conductivity  of  the  frozen  material  above  the  base  lens 
Kb  =  soil  thermal  conductivity  below  the  lens. 

Equation  20  expresses  the  fact  that  conductive  heat  fluxes  on  opposite  sides  of  the  lens  base  differ 
by  an  amount  proportional  to  the  rate  of  phase  change  as  liquid  flows  into  the  lens.  This  equation 
is  included  with  the  governing  equations  in  eq  14  as  a  simultaneous  condition.  As  heave  proceeds, 

8  is  continually  increased,  both  because  the  location  of  the  freezing  zone  advances  towards  a  stready 
state  and  because  heaving  transfers  material  to  the  zone  above  the  lowest  lens. 

Another  boundary  condition  at  x  =  b  is  provided  by  a  force  balance  there.  On  the  lens  side  of 
the  boundary,  the  ice  is  subjected  to  the  full  overburden  pressure  P.  Although  ice  pressure  below 
the  lens  varies  from  P,  it  is  nevertheless  continuous  across  the  lens  boundary,  so  that  u-Jx-b)  is 
still  equal  to  P.  Combining  this  information  with  eqs  1  and  2  provides  an  equation  that  serves  as  the 
necessary  second  condition  at  x  =  b: 

P-u  -  =  toiw(/4u  +  BT)  .  (21) 

This  relation  is  also  incorporated  into  the  governing  set  of  equations. 

An  important  quantity  at  this  stage  in  the  analysis  is  the  neutral  stress  on,  which  is  taken  to  be 
a  weighted  combination  of  isotropic  stresses  in  the  two  phases  of  pore  constituent: 

o„  =  XU+  (1  -  x)«j  •  (22) 

The  quantity  x(0,}  is  an  empirically  determined  stress  partition  parameter,  specifying  the  relative 
participations  of  u  and  in  the  resultant  on.  As  freezing  progresses  below  x  =  b,  ice  content,  \jj 
and  </j  increase.  At  the  same  time,  the  value  of  u  near  the  lens  is  held  down  in  accordance  with 
eq  21.  In  general,  on  passes  through  a  maximum  somewhere  beneath  x-b  and  eventually  surpasses 
the  value  of  P.  When  this  happens,  the  pore  contents  alone  are  able  to  support  the  overburden, 
no  supportive  force  is  transmitted  through  the  soil  grains,  and  a  new  lens  forms.  In  practice,  on 
was  considered  to  surpass  P  when  it  reached  a  critical  value  otr,  approximately  1%  greater  than  P. 

With  the  formation  of  a  new  lens  at  the  location  of  ocr,  the  numerical  mesh  is  shifted  downwards, 
so  that  x  =  b  coincides  with  the  bottom  of  the  new  lowest  lens.  The  above-freezing  boundary  re¬ 
mains  fixed.  The  analysis  is  repeated,  eventually  another  lens  forms,  the  mesh  is  shifted  again,  and 
so  on. 

The  third  and  last  stage  begins  when  Ts  no  longer  decreases,  in  this  case  reaching  a  minimum  of 
-0.5°C.  The  analysis  proceeds  in  the  same  manner,  as  all  effects  slow  down,  asymptotically  approach¬ 
ing  a  steady  state.  Eventually  no  new  lenses  form,  and  the  lens  growth  rate  V\  approaches  zero. 
During  this  and  other  stages,  the  cumulative  heave  h  at  any  time  t  is  V^dt,  which  may  be  evalu¬ 
ated  progressively  through  time  by  simple  quadrature. 


CALCULATED  RESULTS 

As  a  test  of  reasonableness,  simulations  were  performed  over  all  three  stages  for  a  hypothetical 
soil  column,  as  described  above.  The  values  for  2(pc0)n  and  Kh  were  assumed  to  be  2.0  x  10® 
J/°C  m3  and  4.0  J/m  s  °C,  respectively,  where  ice  was  present,  and  2.8  x  106  J/°C  m3  and  3.0 
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j/m  s  °C  where  it  was  not.  Gravity  was  neglected.  (Neither  the  theory  nor  the  computer  program 
are  invalidated  if  these  assumptions  are  discarded.)  The  hydraulic  conductivity  £(0t)  was  assumed 
to  be  close  in  value  to  that  in  an  unsaturated  soil  with  the  same  liquid  water  content.  It  was  rep¬ 
resented  as  (cf.  Bresler  et  al.  1 978) 

•w-wC-Sr*)’ 

where  the  porosity  0O  is  assumed  to  be  0.4.  A  convenient  relation  for  Q^)  was  adopted  (Fig.  1 ), 
suggested  by  the  data  presented  by  Horiguchi  and  Miller  (1980).  The  experimental  results  obtained 
by  Snyder  (1980)  led  us  to  assume  that  *(0,)  could  be  represented  plausibly  by  (1  -8  J80)'-S. 

Figure  2  shows  dT^,  the  time  elapsed  between  successive  lens  initiations,  as  a  function  of  time. 
At  t3,  Ts  reaches  -0.5°C  and  the  third  stage  begins,  during  which  dT L  increases  markedly.  The 
point  furthest  out  in  time  represents  the  last  lens  initiation  simulated,  at  a  time  when  the  zero  de¬ 
gree  isotherm  had  nearly  ceased  penetration.  Calculated  h  versus  time  is  shown  in  Figure  3,  to¬ 
gether  with  the  ice  ratio  rt,  for  times  when  new  lenses  were  initiated.  Here  rt  equals  h  divided  by 
the  current  length  of  the  frozen  zone  8.  The  magnitude  of  h  increases  throughout,  and  both  quan¬ 
tities  increase  during  the  third  stage.  During  that  stage  each  quantity  increases  at  a  slower  rate  on 
an  arithmetic  time  scale. 

Similar  quantities  may  be  determined  for  each  layer  in  the  system.  With  the  inception  of  a  new 
lens,  another  layer  is  added  to  the  relatively  inactive  zone  above.  That  layer  consists  of  a  section 
of  frozen  soil  above  the  new  lens,  together  with  the  lens  that  had  just  grown  above  that  soil.  The 
final  thickness  of  the  old  lens  hL  divided  by  the  total  layer  thickness  provides  a  layer  ice  ratio  zjL. 
Figure  4  shows  thatr|L  and  /»L  remain  relatively  constant  during  the  second  stage,  increasing  by 
about  one  (rjL)  and  three  (/»L)  orders  of  magnitude  through  stage  three. 

Typical  profiles  of  u,  u{,  on  and  T at  a  point  in  time  in  stage  three  are  shown  in  Figure  5.  At 
this  time,  the  maximum  on  has  just  surpassed  ocr,  and  a  new  lens  will  form  where  that  maximum 


Figure  I.  Ice  content  as  a  function  of  \li  and  as 
a  function  of  T  when  u  =  0. 
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Figure  2.  Time  elapsed  between  successive  lens  ini¬ 
tiations  dTt  as  a  function  of  time. 


Time  (») 


Figure  3.  Calculated  cumulative  heave  h  and  ice  ratio  r-  for 
times  when  new  lenses  were  initiated. 


exists.  The  liquid  pressure  u  drops  sharply  across  the  freezing  zone.  During  earlier  stages  freezing 
is  more  rapid  as  the  surface  temperature  is  lowered,  and  the  relatively  rapid  formation  of  ice  causes 
the  liquid  pressure  to  build  up.  Gradually  this  pressure  is  relieved  on  the  freezing  side  as  liquid 
flows  towards  the  surface  in  conjunction  with  heaving  and  on  the  unfrozen  side  as  water  is  ex¬ 
pelled  from  the  column.  This  leaves  a  maximum  in  between,  which  eventually  subsides  into  the 
kind  of  configuration  illustrated  in  Figure  5. 
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Distance  Below  Lowest  lens  (cm) 


a.  Lens  thickness  ht  and  layer  ice  ratio  rjL  vs  time. 


b.  Lens  locations  (shaded)  and  thicknesses  are  shown 
schematically  in  the  portion  of  a  column.  The  numbers 
next  to  the  lenses  are  their  initiation  times. 


Figure  4.  Calculated  distributions  of  segregated  ice. 


Pressure  (hPo) 


Temperature  ( *C) 


Figure  5.  Profiles  of  u,  u,  o„  ond  T  at  1.38 x  10s  s. 
The  neutral  stress  o„  has  Just  surpassed  the  critical 
value  o,_ 


SOME  PERSPECTIVE  ON  THE  MODEL 


It  is  gratifying  to  see  simulations  emerge  from  the  mathematical  model  which  are  self-consistent 
and  which  are  qualitatively  consistent  with  common  sense,  with  alternative  calculations  for  simpler 
cases  (Miller  and  Koslow  1980),  and  with  general  experience  with  real  heaving  systems.  Yet  ques¬ 
tions  remain.  The  model  is  based  on  the  idea  that  thermally  induced  regelation  imposes  require¬ 
ments  for  hydraulic  flow  in  the  frozen  fringe  and  that  these  requirements  are  adequately  achieved 
in  interactions  between  macroscopic  equations  for  hydraulic  conductivity  and  thermal  conductivity. 
The  details  of  the  microscopic  recirculation  of  both  water  and  heat  associated  with  regelation  are 
thereby  left  moot  (cf.  Miller  et  al.  1 975).  Although  the  approach  here  provides  a  sufficient  set  of 
equations  to  obtain  a  solution,  it  may  leave  out  interactions  that  should  be  included.  In  other 
words,  it  is  not  altogether  clear  to  us  whether  such  microscopically  based  corrections  are  implicit 
in  the  macroscopic  equations  or  whether  virtual  counterflows  of  water  and  heat  must  be  super¬ 
imposed  as  corrections  to  the  macroscopic  equations.  If  significant  recirculations  have  not  been 
included  adequately,  the  real  “heaving  engine"  may  not  function  as  effectively  at  high  rates  of 
heave  as  the  model  assumes.  Perhaps  these  considerations  will  prove  unnecessary,  but  it  would  not 
be  altogether  surprising  if  there  were  systematic  discrepancies  between  simulations  that  bypass 
these  concerns  and  actual  heave  test  data. 

Meanwhile,  efforts  to  resolve  the  questions  that  lurk  in  our  minds  go  forward.  A  recent  contri¬ 
bution  by  J.R.  Philip  (1980)  should  be  of  great  value  in  estimating  relationships  between  micro¬ 
scopic  and  macroscopic  thermal  fields  associated  with  regelation.  Before  that  contribution  can  be 
exploited,  however,  it  will  be  necessary  to  carry  out  a  companion  analysis  of  hydraulic  impedance 
to  regelation  (a  factor  deliberately  excluded  from  Philip's  analysis)  and  to  assemble  these  and  the 
macroscopic  function  into  an  integrated  model.  Meanwhile,  it  seems  urgent  lo  establish  a  series  of 
experiments  that  are  controlled  in  a  manner  that  can  be  readily  simulated  with  the  computer  pro¬ 
gram  as  it  stands  or  with  modifications  mandated  by  further  thought  and  experience. 
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